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Abstract 

An absorbing non-reflecting boundary condition (NRBC) for practical computations in fluid dynamics 
and aeroacoustics is presented with theoretical proof. This paper is a continuation and improvement of 
a previous paper [2] by the author. The absorbing NRBC technique is based on a first principle of non- 
reflecting, which contains the essential physics that a plane wave solution [1] of the Euler equations 
remains intact across the boundary. 

The technique is theoretically shown to work for a large class of finite volume approaches. When 
combined with the hyperbolic conservation laws, the NRBC is simple, robust and truly multi-dimensional; 
no additional implementation is needed except the prescribed physical boundary conditions. Several 
numerical examples in multi-dimensional spaces using two different finite volume schemes are illustrated 
to demonstrate its robustness in practical computations. Limitations and remedies of the technique are 
also discussed. 

1 Introduction 

It is well-known that non-reflecting boundary conditions (NRBCs) play a key role in fluid dynamics and aeroa- 
coustics computations and still remain an important topic in the current research areas. A spurious reflection resulting 
from an inappropriate numerical boundary condition contaminates the flow field and may eventually spoil the entire 
flow computation. Research on NRBC is a challenging topic in engineering and applied mathematics. A vast number 
of papers on the topic of NRBCs have been published in past decades to treat various computational problems. Details 
can be found in the recent reviews [3, 4, 5] and the references cited therein. 

Generally there are two types of non-reflecting boundaries (see e.g., [3]). One (Type I) is the true domain boundary 
where physical boundary conditions are specified and the other (Type II) is the artificially imposed boundary. 

In this paper, we focus only on the Type I absorbing non-reflecting boundary. For a true boundary where the 
actual physical boundary conditions are given as per the physical problem, the NRBC needs to play a dual role: on 
the one hand it must be consistent with the given physical boundary conditions, on the other hand, it must be able to 
absorb waves or disturbances from the domain interior. 

For the absorbing NRBC, there are two ways to absorb the outgoing waves/disturbances; the ‘hard absorbing’ and 
the ‘soft absorbing’. In the ‘hard absorbing’, the given boundary condition must be modified to an admissible one. 
This can be done using the characteristics method. In one-dimensional flow, at the boundary, Enquist and Majda [6], 
and Hedstrom [7] showed that a boundary condition considered to be non-reflecting is equivalent to saying that the 
characteristic variables corresponding to the incoming characteristic curve remain unchanged across the boundary (see 
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also Hirsch [ 8 ], p. 370). When the flow is subsonic, Riemann invariants (one-dimensional characteristic variables) 
are used to turn a prescribed boundary condition into an admissible NRBC. For multi-dimensional flow, this one- 
dimensional technique is combined with dimension splitting. Such a combined treatment has been the topic of many 
papers (see e.g. [9] and references in [3]). The characteristics based NRBC is the most popular numerical treatment, 
despite its complicated implementation. 

The ‘soft absorbing’ treatments found in the literature generally include (but not limited to): 

(i) A buffer zone (sponge zone) technique can be used to add a non-physical zone between the domain core and the 
boundary to diminish the strength of the outgoing waves/disturbances before they reach the boundary, thus minimizing 
the reflecting effect. In the buffer zone, either artificially increased dissipation is applied or grid stretching is imposed 
(see e.g. [ 10 , 18]), which is equivalent to numerical filtering. 

(ii) In the recently developed perfectly matched layer (PML) method [11, 12] , a specially designed equation 
system is imposed in the matching layer (or absorbing layer) at the boundary to guarantee the exponential decay of the 
disturbances in the layer. 

In addition to the existing absorbing NRBCs, as mentioned in [2], a new absorbing (Type I) NRBC technique was 
recently discovered in an empirical way. The technique is robust but simple to implement and has been successfully 
employed in many computations [16, 17, 18]. 

It is the purpose of the present paper to revisit the basic theoretical background of the absorbing boundary condi- 
tion and to explain why the technique works well for flows in one to three dimensional spaces. Starting from a first 
principle of non-reflection, i.e., there is no reflection at any interior point of a computational domain, a C 1 continuity 
criterion for flow variables is introduced as a generic NRBC. The criterion is then incorporated in the time-dependent 
hyperbolic conservation laws in discretized form, and the consequent absorbing (Type I) NRBC technique is described 
and theoretically proven. 

In the present paper, we first develop the basic principle and criterion for the continuous problem and then apply 
them to the discretized numerical problems. The paper is arranged as follows: in §2, based on a first principle of 
non-reflection, the C 1 continuity criterion of NRBC is introduced and proven. Based on this continuity criterion, the 
absorbing NRBC is depicted in §3. Various numerical examples in two or three-dimensional spaces with two different 
finite volume schemes are illustrated to demonstrate the capability of the new NRBC. Limitations and remedies of the 
absorbing NRBC are discussed in §5. Finally, the paper is concluded with remarks in § 6 . 

Throughout the present paper, 

V = (v 1 ,v 2 ,v 3 ,v 4 ,v 5 ) t = (p,u,v,w,p) T 


is employed to represent the primitive flow variables (for notations, see below). The primitive flow variables V and 
the original non-conservative form of the Euler equations are conveniently employed in the theoretical investigation 
but hyperbolic conservation laws in discretized form are needed in the numerical implementation. Here is the standard 
three-dimensional conservation form of the non-dimensional Euler equations: 

Ut + F x + G y + H z = Q, (1) 


where x, y, z and t are the streamwise and transverse coordinates and time, respectively. The conservative flow 
variable vector, U, and the flux vectors, F, G, H, are given by: 


U = 




( p \ 


( f ^ 


( PV \ 


/ Pw \ 

U 2 


pu 


pu + p 


puv 


pwu 

u 3 

= 

pv 

,F = 

puv 

, G = 

pv 2 + p 

,H = 

pwv 

u A 


i pW i 


puw 


pvw 


pw 2 + p 

\uj 


\ pe / 


\ puH / 


\ pvH ) 


\ pwH / 


where u,v,w and p. p are respectively the velocity components, density and pressure, with the internal energy e = 
+ 1/2 (it 2 + v 2 + w 2 ), the enthalpy H = p/p + e, and 7 = 1.4. The right hand side Q is the source term which 
may possibly include the external forcing terms. 

By considering (x, y, z, t) as coordinates of a four-dimensional Euclidean space, E 4 , and using the Gauss diver- 
gence theorem, it follows that Eq. (1) is equivalent to the following integral form of the conservation laws: 



QmdV", 


to = 1, 2, 3, 4, 5, 


( 2 ) 


NAS A/TM— 2005-2 1 3654 


2 



where S(V) denotes the surface around a control volume, V, in £4 and I m = (F m , G m , H m , U m ) represents the flux 
vectors. The vector differential ds = nds, where n is the outgoing unit normal vector of the infinitesimal surface 
element ds in £4. 

2 A first principle of non-reflection and the continuity criterion 

In this section, a theoretical analysis of non-reflecting is conducted for the continuous problem (as opposed to a 
discretized problem) of the Euler equations. Then, it is applied to the discretized NRBCs in the sections that follow. For 
clarity, the flow variables are assumed to be continuous near the boundary, and all the spatial and temporal derivatives 
exist. It should be stated that, even though discontinuities such as shocks or contacts may develop and propagate across 
the boundary, for a discretized absorbing NRBC, such continuous flow requirement is always satisfied locally at the 
boundary via the numerical procedure, and hence the following first principle and criterion of non-reflecting are still 
valid. 

From the classical theory (e.g. Courant and Hilbert[l]) of partial differential equation (PDE), it is well known that 
hyperbolic PDE systems support wave solutions. In the present paper, we specifically focus on the Euler equations. 
Due to their hyperbolicity, they support linear continuous waves such as pressure pulse waves, entropy (density) waves, 
and vorticity waves as well as nonlinear or discontinuous waves such as shocks. 

2.1 A first principle of non-reflection 

As the PDE system supports wave solutions, the next issue is focussed on the questions: what is ‘non-reflecting’ 
and how a boundary can become non-reflecting? Hereby, instead of following the conventional wisdom of a sig- 
nal/wave propagating along characteristics, the concept of ‘non-reflecting’ is viewed from a different but simple per- 
spective. 

Recall a general observation that, no matter how crude the numerical scheme is, spurious reflection can only 
be generated at the domain boundaries, no reflection ever occurs at any interior point. This is because only at the 
boundary is the numerical solution subjected to two competing conditions, namely, the governing equations and the 
imposed boundary conditons. 

The following logical argument further justifies that no reflection should occur in the domain interior: At an 
interior point of the domain, the solution is subjeted only to the governing equations. For any scheme that converges, 
the numerical result at this point, which may even contain the possible spurious reflection propagated from a boundary, 
is literally considered to be 100% a discretized solution of the governing equations and nothing else, disregarding any 
discretization error. Conceptually, no such ‘impurity’ as ‘reflection’ is generated at this interior point. 

Theoretically, for a hyperbolic PDE system, locally in a small neighborhood of any interior point O, by using 
Cauchy’s method of the Fourier integral (e.g.,[l], p.210, or [ 8 ]), the homogeneous part of the solution of the initial 
problem, which may contain oscillations, can be expressed by the superposition of plane waves: 

M />oo roo 

■■■ V m (k )e ie dk 1 dk 2 ...dk n , (3) 

m=1 J-oo J - 00 

Here, M is the number of PDE equations in the system, k = (k \ , ky . ..., ) is the /(-dimensional wave number vector, 

with kj, j = 1, 2, ..., n, ranging from —00 to 00, i = \J — 1, and x = (xi,X 2 , ■■■, x n ) is the ;r-dimensional position 
vector of this point O. V m = V m (k) and 6 = 0(k, x, t) = x • k — Lu m t are respectively the amplitude and the phase 
of the plane wave, where uj m = o> m (k) (a real number) is the angular frequency and can be determined as a function 
of the wave number vector k by substituting the single plane wave Ve ie into the PDE system. For each plane wave 
component in (3), as its amplitude and phase are functions of k and remain intact when passing through this interior 
point O, no distortion or reflection occurs. So also for the solution V. 

Based on such observation, logical reasoning, and the theory of hyperbolic PDE system, we conclude that if 
all the points at the boundary are actually interior points of the domain, there is no reflection. A first principle of 
non-reflection is then proposed and adopted in the present paper: 

A boundary of the solution domain is said to be non-reflecting if every point on the boundary surface can be 
considered as an interior point. 

There are two requirements defining an interior point: 
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Figure 1 . a: Sketch of a one-dimensional scalar plane wave, b: the curves v = v(x, t) on a plane t = to at the boundary x = xq. showing 
intuitively that C 1 continuity is necessary for non-reflection. 


1. the point must be physically located within the domain; 

2. the consistency requirement that the governing PDE system is the only condition to be satisfied at this point. 

Other conditions, if any, must be equivalent to or consistent with the PDE system. 

In order that all the points along the boundary be interior points, of course, the solution domain must be extended 
to include at least a small neighborhood around the boundary. Numerically, the ghost cell is a natural choice for the 
neighborhood at the boundary. The consistency requirement is a more important issue. An NRBC is sought such that 
it fully complies with the governing equations, that is, solutions in the domain interior and in the domain exterior must 
be mutually a solution continuation of each other, satisfying the same PDE system in the extended domain. 

This concept or definition of non-reflection leads naturally to a generic NRBC - the continuity criterion for Euler 
equations, which will be proposed and proven in §2.2. The definition is a general concept and may be applied to 
hyperbolic systems other than the Euler equations. 


2.2 The continuity criterion - a generic NRBC 

With the presence of the boundary ( a cylindrical hyper-surface), the four-dimensional Euclidean time-space E 4 is 
bisected into two portions, domain interior Di and domain exterior D e (Figure la illustrates a 1-D scalar plane wave in 
E 2 , only the real component of V is shown. The subscripts i and e stand respectively for domain interior and domain 
exterior ). Within Di, the flow is governed by the original Euler equation system. It is also assumed that in the domain 
exterior D e , at least in a neighborhood of the boundary, the flow is still governed by an Euler equation system. If the 
two Euler equation systems are identical at all the points of the boundary, the consistency requirement is satisfied, and 
these points are interior points of the extended domain. The boundary is hence non-reflective. 

The continuity criterion leads naturally to such non-reflection, it is stated as follows: 

At each point of the boundary, the flow variables V and their spatial derivatives, namely, ~V X , V y , V 2 , are 
continuous; or, V is C 1 continuous in space across the boundary. 

It can be easily shown that locally there is no reflection at the boundary surface if the continuity criterion is 
satisfied. 

Proof: Consider the non-conservation form of the Euler equations ( e.g., see [8]) for both domain interior and 
exterior: 


dV ~<9V 

dt dx 


B 


dV 

dy 


~dV 

C— = Q. 

oz 


(4) 


Here Q is a nonhomogeneous source term vector and may be a vector function of x, y, z, t and V, and could include 
forcing terms. 

Note that the jacobian matrices A, B, and C are directly functions of V ( details may be found in [8]). As a result 
of the C 1 spatial continuity of V, all the terms (items) 


dV dV dV 

~dt ’ Hx' Hy' 


dV 

Hz' 


A,B,C , and Q 


in the governing Euler equations are identical across the boundary surface. Continuity of is directly inferred 
from (4). Now, (4) can be considered as a single governing equation system in an extended domain which covers 
the boundary as its interior. The consistency requirement is satisfied at all the boundary points; hence, there is no 
reflection at the boundary. 
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2.3 Physical interpretation of the continuity criterion 

For the three-dimensional Euler equations, corresponding to Eqn. (3), M = 5 and n = 3. As a special case of the 
hyperbolic PDE system in §2.1, locally in a small neighborhood of an interior point, the homogeneous solution of the 
Euler equations can always be decomposed into plane waves having a vector form [8]: 

V = Ve ie , (5) 

Here, as in §2.1, V = V (k, x , y, z) is the ‘amplitude’, and 6 = 0(x, t) = x • k — uit is the phase of the plane wave, 
where k = (k x , k y ,k z ) is the given wave number vector and to = w(k) (a real number) is the angular frequency, and 
x = (. x , y , z) the position vector of the point O. Note that 0 is a real number; 9 = const, represents a space-time 
characteristic surface or wave front [1], 

Now, as the continuity criterion guarantees that a boundary point can be considered as an interior point, the plane 
waves are employed in this paper to serve two purposes: 

1. to demonstrate the physical effect of the continuity criterion: at such a boundary point, a plane wave, i.e., its 
amplitude and phase remains intact or undistorted across the boundary. Hence the non-reflecting property is 
justified as stated in §2.1; 

2. to conduct phase error analysis at the boundary (in §5.1). 

Figure lb sketches how a simple one-dimensional scalar plane wave propagates across a boundary x = Xq in a 
plane t = to. The upper curve corresponds to C 1 continuity and exhibits no reflection. The lower curve corresponds 
to C° continuity with a possible kink of slope or phase error, spurious reflection may occur. 

The C 1 continuity criterion is a generic NRBC. When incorporated in certain finite volume approaches with the 
hyperbolic conservation laws, a new class of absorbing NRBCs is created. As will be shown in §3 and §4, these 
NRBCs are simple but robust. The continuity criterion may also have direct applications, but that would be the topic 
of another NRBC paper. 

2.4 Relation of the continuity criterion to the characteristics-based NRBCs 

The continuity criterion is closely related to the well-known characteristics based NRBCs [6 - 9], although seem- 
ingly there is no apparent way to infer one from the other. First, both rely on the hyperbolicity of the equation system, 
which supports wave/disturbance propagation and characteristics. Secondly, for the characteristics-based NRBCs, as 
suggested in [6, 7], for the incoming characteristic variable Wk, Awjfc should be set to zero across the domain boundary. 
This step is consistent with the continuity criterion too. As such, their non-reflecting performance is expected to be 
similar. 

However, the continuity criterion treats a multi-dimensional plane wave as a physical entity and keeps it intact or 
undistorted across the boundary surface to attain the ‘non-reflecting’ effect. This contrasts to the usual treatment of 
characteristic field decomposition in the characteristics-based NRBCs and is much simpler. 

3 The absorbing NRBC and the absorbing layer (matched layer) 

The continuity criterion is primarily applied to the absorbing NRBC. This absorbing NRBC is generally con- 
sidered to be a difficult problem because the boundary conditions are subject to excessive restrictions. On the one 
hand, the given flow conditions must be satisfied; on the other hand, the boundary must be non-reflecting or absorbing 
for the outgoing waves. Conventionally, by determining the one-dimensional outgoing or incoming characteristics, 
the characteristics-based NRBCs modify the boundary conditions to an admissible one that handles the absorbing 
boundary problem. The process is termed ‘hard absorbing’ in §1 and becomes cumbersome when the flow is multi- 
dimensional and complicated. The perfectly matched layer or PML method [11, 12] offers an attractive treatment, 
in which a set of equations in the matched layer is prescribed and solved such that the outgoing wave is dissipated 
exponentially and absorbed. Such absorbing occurs in the matched layer and hence belongs to the ‘soft absorbing’. 

The new technique is shown to work for a large class of finite volume schemes that follow the two-stage procedure: 
an evolution (time-marching) stage and a reconstruction stage. The continuity criterion is here combined with the 
conservation laws in discretized form and yields a desirable absorbing NRBC - simple but robust. The non-reflecting 
property of the Type I absorbing NRBC is established in §3.1. Section 3.2 provides a general guideline for numerical 
implementation of the absorbing NRBC. The matched layer or absorbing layer is discussed in §3.3. 
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boundary surface element AS 




Figure 2. a: Boundary surface element AS as an internal interface of the control volume V ; b and c: Examples of2-D numerical implementation: 
b: A standard upwind scheme with an unstructured triangular grid; ‘R ’ and ‘U represent respectively the right and left states for the Riemann solvers, 
they are obtained by extrapolated flow data from the current cell center O and the neighboring cell centers D, E and F. The prism cylinder based 
on A ABC is the space-time control volume. The boundary surface element AS practically coincides with the side AB. c: For the CE/SE method, 
with AB as the boundary AS, the space-time control volume V is a hexagon cylinder based on ADBECF A. In both cases, physical boundary 
conditions are given at the ghost cell center D. 


3.1 The absorbing NRBC 

After discretization, the boundary surface S is also divided into small boundary elements, A,S', according to the 
grid. Let V be a control volume which contains the boundary surface element AS as its internal interface (Fig. 2a). 
Note that V may or may not be a grid cell. It will be shown that AS is non-reflective with any prescribed physical 
boundary conditions at the ghost cell center. No additional numerical implementation is needed. The mechanism 
works as follows: the C 1 continuity of the flow variables V across AS provides the non-reflecting property as shown 
in the next paragraph. Then the inconsistent flow data between the imposed boundary conditions and the outgoing 
wave/disturbance is treated as a ‘numerical jump’ around the boundary, captured, and absorbed by the shock-capturing 
scheme which is based on the hyperbolic conservation laws. As a result, a matched layer (ML) or absorbing layer is 
automatically generated around the boundary. 

It is straightforward to prove the ‘non-reflecting’ property of the boundary surface element AS. In a finite volume 
numerical approach, at the evolution stage (time-marching stage), for any prescribed physical boundary conditions at 
the ghost cell center, the integral conservation laws (2) are solved with the control volume V, yielding a spatial average 
V (or U) over V at the new time level n + 1 . At the reconstruction stage, by using a flux limiter, such as the minmod 
or the van Albada limiters etc., to recover the spatial gradients of V (i.e., V x , V y and V 2 ), this constant averaged V 
is rebuilt into V r,+1 - a linear, high resolution (or second order accurate ) solution of (2) within the control volume 
V. It is the reconstruction stage that provides the non-reflecting property. Note that V rl+ I represents a discretized 
solution of (2) or a discretized ‘weak solution’ of (4) within the entire control volume V and that AS is an internal 
interface of V. As a linear vector function within V, V” +1 satisfies the C 1 continuity criterion at the surface center 
of AS (Fig. 2a), hence AS is non-reflective. 

3.2 Implementation of the absorbing boundary conditions 

Construction of the absorbing NRBC is straightforward too - no additional numerical implementation is required 
except the prescribed boundary conditions at the ghost cell centers, because an absorbing layer is automatically gen- 
erated. For the class of finite volume approaches which contain the two-stage procedure of evolution (time-marching) 
and reconstruction, here are some general guidelines: 

1 . Choose an appropriate control volume, V, at the boundary and ensure that the boundary surface element. A, S', is 
an internal interface of V. For most of the finite volume schemes, usually, a control volume V where the Gauss 
divergence theorem (2) is applied is equivalent to a grid cell. Then, AS may coincide with the cell surface at the 
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Figure 3. Comparison of absorbing NRBCs at an inlet boundary: a: Conventional ‘hard’ absorbing NRBC (characteristics-based); b: Type I ‘soft’ 
absorbing NRBC, the matched layer plays the role as a buffer layer between the inconsistent flow data at domain boundary and interior. 


domain boundary but is considered lying at the inner side of the boundary surface (Fig. 2b). For some schemes, 
e.g. the CE/SE scheme [14,15], a control volume V is larger than a grid cell and AS is truly an internal interface 
of V (Fig. 2c). 

2. The finite volume approach follows the standard two-stage procedure: the evolution (time-marching) stage and 
the reconstruction stage by flux limiter. Nevertheless, even when the reconstruction stage is missing, the NRBC 
still holds, as the case is equivalent to a first order scheme, with all the spatial gradients vanished. 

This absorbing NRBC is flexible and valid with practically any cell shapes, structured or unstructured. Figure 2b 
and 2c show examples of the two-dimensional numerical implementation - the treatments in a standard upwind scheme 
and in the CE/SE method. Both will be used as numerical platforms for the computed examples. Flow data at the 
boundary cell center O is updated based on flow data at its three neighboring cell centers D, E and F. The divergence 
theorem is applied to the control volumes: the triangle cylinder A ABC or the hexagon cylinder ADBECF. The side 
AB is the non-reflecting boundary and an internal interface of the control volume as well. 

3.3 The matched layer (absorbing layer) near the boundary 

As the specified boundary conditions are imposed at the ghost cell centers, they also play a dual role of an 
absorbing (Type I) NRBC. Often, after a considerable long run-time, inconsistency between the imposed boundary 
conditions and the flow in the domain interior near the boundary may develop. Then, while the C 1 continuity across the 
boundary surface rules out spurious reflections, the shock-capturing scheme automatically captures the inconsistency 
in a ‘matched layer’ (ML) or ‘absorbing layer’ . Decaying of the inconsistency in the absorbing layer is fast. Depending 
on the required accuracy, the matched layer (ML) is typically only a few cells thick. It is non-physical, but plays an 
important role in separating and matching the inconsistent flow data at the boundary and the domain interior. The 
matching or absorbing is thus considered ‘soft’. 

The situation is similar to that in the PML (perfectly matched layer) method [11, 12], except that the ‘matched 
layer’ or absorbing layer is automatically generated and no special equations in the matched layer are needed. 

4 Numerical examples 

The absorbing NRBC is the primary but indirect application of the continuity criterion. As a matched layer or 
absorbing layer will be formed automatically, the absorbing NRBC is flexible and capable of handling complicated 
non-reflecting situations, such as with the ‘difficult’ cases where locally the flow disturbance may enter or exit the 
domain through the boundary and where the flow conditions could be supersonic along a portion of the boundary but 
subsonic along another portion of the same boundary. 

Since computations in multi-dimensional space are of primary importance, in this section, the absorbing NRBC 
is tested in examples in two- or three-dimensional space with different finite volume schemes to demonstrate its 
performance and robustness. Here, a standard upwind scheme with an unstructured triangular grid, as depicted in 
Fig. 2b (for 2-D problems) and the space-time conservation element and solution element (CE/SE) method [13, 14] 
are chosen as the numerical platforms for computing the examples. The upwind scheme uses Roe’s approximate 
Riemann solver and renders first order accuracy in time but high-resolution in space. The CE/SE scheme is nominally 
second order accurate in both space and time. 


NAS A/TM— 2005-2 1 3654 


7 


vorticity/ 

entropy 

wave 



t=30 





Figure 4. Propagation of a linear acoustic pulse with vorticity/entropy disturbances: density contours plots att = 6, 30, 60; with comparison to 
exact solution along the x axis, showing no reflection at the outflow boundary. The CE/SE scheme is employed here. 




Figure 5. Plots 1-12, propagation of a Gaussian pulse in a circular domain from t = 60 to t = 120. The flow Mach number, M, is 0.2 at an 
angle of 45°. Although the flow is inclined and the boundary is curved, no visible reflection is observed. A first order Godunov upwind scheme is 
employed in this problem. The absorbing layer exists but probably too thin and weak to be visible. 


4.1 2-D acoustic pulse, entropy wave, and vorticity wave propagation 

This problem illustrates the propagation of the three basic types of linear waves in a two-dimensional domain [15, 
16]. The computational domain in the x-y plane is a square with —100 ^ x ^ 100, and —100 ^ y ^ 100. A uniform 
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Figure 6. 




Plots 1-6, propagation of a Gaussian pulse in an egg-like domain at t = 150 , 200 , 250 , 300 , 350 , and 400 . The flow Mach number, 


M, is 0.05 and directs at 45 °. ( about 65,000 unstructured tiiangular cells, with the CE/SE scheme) 


grid with 40, 000 uniform triangulated cells is used with the CE/SE scheme. Initially, a Gaussian acoustic pulse is 
located at the center of the domain (0,0) and a weaker entropy/vorticity disturbance is located off the center at (67,0), 
with a mean flow Mach number, M = 0.5, in the x direction. The corresponding non-dimensional pressure, density, 
and streamwise and transverse velocity components are given by 

p=- + Se~ a ^ x2+y2 \ p= 1 + Se~ aiix2+v2) + 0.1 Se- a2 ^ x - 67 ^ +y2] , 

7 

u = M + 0MSye- a2[( ' x - 67)2+y2 \ v = -0.045(a; - 67)e~ a2[{x - 67)2+y2 \ 

where a\ = In 2/9, 0.2 = In 2/25, and with a small amplitude factor, S = 0.001, the two-dimensional Euler equations 
are practically linearized. At all four boundaries, flow variables are given with the M = 0.5 mean flow conditions: 

p a = — , p = l,u = M = 0.5, v = 0. 

7 

These conditions play the role as an absorbing NRBC. Figure 4 illustrates the evolution of the density contours at 
different time steps with comparison to the exact solutions (perturbation part) along the x axis. It is observed that the 
numerical solutions agree well with the exact solutions and that there is no reflection along x axis. As the the waves 
are weak, the absorbing layer is not visible. 

4.2 2-D acoustic pulse propagation 

The second example is a two-dimensional subsonic non-linear wave propagation problem. A circular computa- 
tional domain in the x-y plane is used, ranging from —100 ^ x ^ 100, and —100 ^ y ^ 100 (Fig. 5). There are about 
87,200 unstructured triangular cells in the domain, and the two-dimensional unsteady Euler equation system is solved 
by a simple first order accurate Godunov upwind scheme. Initially, a Gaussian acoustic pulse is located at (10, 10) in 
the domain, with a mean flow Mach number, M, of 0.2 in a direction 45° inclined to the x axis. 

p = e ~‘ !n2 ( x2 +V 2 )/9 + 1/1.4, p = g — ^2(z 2 +y 2 )/9 + 1} u = v = 0.1^2. 

The Type I NRBC which is identical to the initial conditions i.e. p = 1/1.4, p = 1, u = v = O.lv/2) is imposed at 
the circular boundary. As time elapses, the non-linear pulse convects with the flow and propagates. Figure 5 shows 
twelve frames of the isobars from t = 65 to f = 120 at an increment of At = 5. When the pulse propagates towards 
the boundary and exits the domain, no visible reflection is observed. The matched layer is too thin to be visualized. 
The curved boundary and the oblique flow direction help to demonstrate the capability of the absorbing NRBC. 
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Figure 7. Two-dimensional radial Riemann problem, showing how a shock and a contact discontinuity are perfectly absorbed by the boundary. 
Plot 1-5 correspond to their propagation at t = 50, 100, 125, 175, and 225. The mean flow is stationary (M = 0). 


4.3 2-D acoustic pulse propagation in an egg-like domain 

This example demonstrates that the absorbing NRBC can be applied to a domain of more general shape. All the 
conditions are similar to the previous example except that the domain has the shape of an egg and that the mean flow 
Mach number, M, is 0.05. The numerical scheme used here is the CE/SE scheme. Figure 6 shows six frames of the 
isobars from t = 150 to t = 400 at an increment of At = 50. In this case, when the pulse propagates towards the 
boundary and exits the domain, the absorbing NRBC performs generally well, except some minor reflection in the last 
3 plots, due to certain phase errors that will be described in §5.1. 

4.4 Propagation of discontinuous waves (1 ) 

Here the problem is a two-dimensional radial Riemann problem and the physical phenomena are complicated. 
This example demonstrates with the CE/SE scheme how a shock and a contact discontinuity are absorbed by the 
boundary. The same circular domain as in §4.2 is employed. The mean flow is stationary: (po,vo,i’o,Po) = 
(0.5, 0, 0, 0.1/1. 4), which are also imposed at the ghost cell centers at the circular boundary as the absorbing NRBC. 
Initially, within an inner circle of radius r = 50, a state of high pressure and density is given: (p,u,v,p) = 
(1.0, 0, 0, 1.0/1. 4). Figure 7 shows five plots of isobars and isopycnics (density contours) at different time steps. 

In Fig. 7a, the first two plots show how the circular shock approaches the absorbing boundary. From the third plot on, 
the shock passes through the boundary and is absorbed by the boundary, leaving no trace of reflection. Figure 7b not 
only demonstrates how the shock leaves the domain, but also the evolution of the contact discontinuity. In the fifth 
plot, the contact already passes through the boundary and no reflection appears. 

4.5 Propagation of discontinuous waves (2) 

This radial Riemann problem is similar to the previous example, except that the mean and core flows move at M = 
0.2 and 45°. Initially, the mean flow and core flow data are respectively (po> uo> vo,Po) = (0.5, 0.1414, 0.1414, 0.1/1. 4), 
and (pi, Ui , Vi, pi) = (1.0, 0.1414, 0.1414, 1.0/1. 4). Figure 8 shows eight plots of isopycnics (density contours) at dif- 
ferent time steps. The circle-like shock moves towards the northeast direction in plots 1 and 2. In plots 3 and 4, the 
shock is almost fully absorbed by the boundary, leaving no visible reflections. From plot 5 on, the contact propagates 
slowly towards the boundary and is well absorbed. A high-resolution upwind scheme with a minmod flux limiter is 
employed for computing this case. 


4.6 2-D subsonic cavity flow 

Figure 9 illustrates the instantaneous isobars for a cavity noise problem. The grid consists of 42, 000 triangulated 
structured cells. The problem is a M = 0.8 flow past a rectangular cavity of aspect ratio of 6.5. At the cavity walls, 
no slip boundary condition is imposed. Due to vortex shedding and acoustic feedback at the cavity edges, strong 
nonlinear acoustic waves are generated and propagate in both upstream and downstream directions [17]. Figure 10 
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Figure 8. Two-dimensional Riemann problem with mean and core flows moving at 45° and M = 0.2, showing how the the shock and contact 
discontinuity are absorbed at the boundary. Plots 1 -8 correspond to the density contours att = 40, 60, 80, 90, 100, 120, 140, and 160. 



Figure 9. Isobars snapshot for a cavity noise problem (Mach number M = 0.8), showing inflow NRBC and its absorbing property. Also, no 
visible reflection is found at the top and the outflow boundary 


is an enlargement of Fig. 9 around the inflow area, revealing the details of the contours at the matched layer. It is 
observed that there is no spurious reflection and that the matched layer (absorbing layer) is about 4-5 cells thick. A 
CE/SE type Navier-Stokes solver is used in the computation. 

4.7 3-D circular jet screech problem 

The absorbing NRBC can be extended to three-dimensional space problems in a consistent way. Figure 1 1 is the 
sketch of an underexpanded circular jet in three-dimensional space. The jet nozzle protrudes into the computational 
domain by S = 2D, D being the inner diameter of the nozzle. The unstructured mesh contains 3.7 million hexahedral 
cells. At the inlet plane, an ambient (stationary) condition or the Type I absorbing NRBC is specified. Jet flow at 
higher pressure (with jet Mach number Mj = 1.42) is specified at the nozzle exit. For over or under expanded jets. 
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Figure 1 0. Details of the contours at the inflow boundary, showing the matched layer (absorbing layer) at the inlet and its spreading over the grid. 
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Figure 1 1 . Left: Sketch of the 3-D circular jet screech problem and computational domain, jet Mach number Mj = 1.42, S = 2 D, L = 10 D 
and R = 5.5 D. Right: instantaneous pressure isosurfaces, showing screech waves, and in particular, the disc of matched layer at the ambient inlet. 


screech tone noise and waves may occur in a feedback loop in which the instability waves (vortices) in the shear 
layer interact with the shock-cell structure in the jet core and generate screech waves that propagate upstream to the 
nozzle lip triggering new vortices. A CE/SE type three-dimensional Navier-Stokes solver is employed for the work. 
Figure 11 also shows a snapshot of the pressure isosurfaces after running 510,000 time steps. The helical screech 
waves are clearly illustrated. At the inlet, an absorbing layer (or ML) is automatically formed as shown in Fig. 11. 

5 Discussions on the NRBCs 

Generally, the new NRBC exhibits little reflection. However, due to discretization errors and the possible conse- 
quent phase error at the physical wave front, there are situations when the NRBC may fail and discernible reflections 
occur. Limitations and remedies are briefly discussed here. 

5.1 Limitation of the NRBCs 

Due to discretization and the possible consequent phase error, the accuracy of the NRBC could be degraded. 
Consider the Fourier mode in a plane wave solution (5), i.e., e 1 ( k#x_a;t ). Here, 0(x, t) = k • x — cot is the phase of 
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Figure 12. Wave propagation across a boundary ; a: wave front coincides with boundary surface, best non-reflecting effect is obtained; b: wave 
front is inclined to boundaiy surface, phase error may occur and cause spurious reflection. Arrows indicate wave propagation direction (direction 
of wave number vector) 



Figure 13. Mach radiation from a M = 2 axisymmetric jet; Top: no buffer zone is added, flow held is contaminated; Bottom: with buffer zone 
added (not shown), showing a clean acoustic held. 


the wave mode, with k being the given wave number vector in the propagation direction, and ui the frequency. The 
characteristic surface, 6*(x, t) = const., stands for a wavefront ( see e.g., Courant and Hilbert [1] or Hirsch [8]). After 
discretization, the center, O, of AS (Fig. 2a) is employed to represent the entire AS (a hyper-surface in both time and 
space), both spatial and temporal phase errors may occur. While the temporal phase error is part of the discretization 
error, the spatial phase error is closely related to the spurious reflections and is less tolerable. Assume that time t is 
held unchanged, the spatial phase error can be easily analysed: Let x = (x, y, z) be the position vector of any point 
on AS and xo the position vector of the center (). Then, the spatial phase error AO due to replacing x by xo is: 

AO = k • x k • xo = k • (x — xo) (6) 

note that the vector Ax = x — xo lies on AS, and AO = 0 when k is normal to AS. Therefore, spurious reflection 
is minimized when the wave front (characteristic surface) coincides with the boundary (Fig. 12a). Otherwise, a phase 
error of order 0( Ax) may be introduced and cause numerical reflection. It should be noted that the above analysis 
involves only discretization and is hence valid for other NRBC techniques such as the characteristics-based NRBC. 

It is conceivable that changing the orientation of the boundary will help to reduce spurious reflections. But in prac- 
tice, an effective remedy is to impose a buffer/sponge zone between the boundary and the domain interior. Although 
the same governing equations (1) or (2) are employed in the sponge zone, numerical damping is highly increased 
to diminish the wave/disturbance amplitude before it reaches the boundary and to minimize the spurious reflection 
typically by grid stretching. An example is illustrated here to demonstrate the effectiveness of the buffer/sponge zone. 
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